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Executive  Summary 


Due  to  a  lower  level  of  funding  than  originally  anticipated,  the  emphasis  of  the  second 
year  of  this  contract  has  been  on  supporting  the  existing  version  of  the  Parameterized  Real-time 
Ionospheric  Specification  Model  (PRISM)  by  providing  updates  and  bug  fixes  as  problems  were 
uncovered  by  Computational  Physics,  Inc.  (CPI)  and  by  other  users  of  PRISM.  Nevertheless, 
work  aimed  at  the  development  of  PRISM  2  has  continued  at  a  low  level.  A  working,  multi-ion 
version  of  GTIM  being  necessary  both  to  this  effort  and  to  a  separate  development  effort  funded 
under  a  separate  contract  (F19628-96-C-0041),  development  of  the  said  version  of  GTIM  took 
place  under  the  other  contract.  However,  the  implications  of  this  work  for  PRISM  2  are 
described  in  this  report  along  with  a  summary  of  the  updates  and  changes  made  to  the  current 
operational  version  of  PRISM. 
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1.  INTRODUCTION 


The  primary  objective  of  this  contract  is  the  improvement  of  the  Parameterized  Real-time 
Ionospheric  Specification  Model  (PRISM).  This  objective  has  several  components  including  (1) 
improvement  of  PRISM’S  base  climatology,  i.e.,  the  Parameterized  Ionospheric  Model  (PIM),  (2) 
extending  PRISM  (and  PIM)  to  include  the  plasmasphere,  and  (3)  improving  PRISM’s  data 
assimilation  algorithm.  The  first  two  objectives  will  be  accomplished  by  parameterizing  a  single, 
global,  physics-based,  multi-ion  ionospheric  model  (AFRL’s  Global  Theoretical  Ionospheric 
Model,  GTIM)  instead  of  the  separate  regional  ionospheric  models  used  in  the  original  PRISM 
development  effort.  The  new  version  of  PRISM  will  be  called  PRISM  2.0  and  will  incorporate 
other  improvements  as  well.  A  secondary  objective  of  this  contract  is  the  continued  support  of 
PRISM  1.6,  which  is  the  operational  version  at  the  50th  Weather  Squadron  at  Falcon  AFB,  (2)  the 
development  of  PRISM  applications  such  as  three  dimensional  ray  tracing  software,  and  (3)  the 
development  of  visualization  software  for  use  with  PRISM  output. 

The  emphasis  of  the  second  year  of  this  contract  was  on  providing  support  for  the  current 
version  of  PRISM  while  preparations  for  the  production  of  PRISM  2  continued  at  a  lower  level 
than  during  the  first  year.  This  was  due  in  part  to  a  difficulty  that  we  encountered  with  the  light 
ion  version  of  GTIM.  This  difficulty  is  described  in  Section  2.  Modifications  to  the  current 
version  of  PRISM  to  correct  problems  that  have  been  uncovered  by  various  PRISM  users  are 
described  in  Appendix  A. 

The  content  of  this  report  reflects  the  state  of  the  effort  as  of  August  1997.  A  number  of 
things  have  changed  since  then,  and  the  reader  is  urged  to  consult  later  reports  in  this  series  for 
more  current  descriptions  of  PRISM  development. 


2.  THE  HIGH  ALTITUDE  EXTENSION  TO  PRISM  AND  PIM 

The  Parameterized  Real-time  Ionospheric  Specification  Model  (PRISM)  has  been 
described  by  Daniell  et  al.  [1990]  and  Anderson  [1993].  The  climatological  model  on  which 
PRISM  is  based,  known  as  PIM,  has  been  described  by  Daniell  et  al.  [1995].  The  validation  of 
an  early  version  of  PRISM  was  described  by  Daniell  et  al.  [1993].  Some  of  the  design 
considerations  for  the  development  of  PRISM  2  were  described  in  Daniell  et  al.  [1998], 

One  of  the  important  innovations  to  be  used  in  this  development  was  a  new  version  of  the 
AFRL  Global  Theoretical  Ionospheric  Model  (GTIM)  that  included  the  light  ions  H+  and  He+, 
allowing  the  modeling  of  the  plasmasphere.  In  the  process  of  implementing  GTIM  on  our  in- 
house  computers  to  make  production  runs  for  PRISM  2  and  for  another  contract,  we  found  the 
code  to  be  subject  to  numerical  instabilities.  In  the  process  of  diagnosing  the  source  of  these 
instabilities,  we  realized  that  there  was  a  problem  with  the  use  of  this  approach  for  the 
plasmaspheric  field  lines  (L  >  ?). 

GTIM  was  developed  from  the  older  code,  LOWLAT,  which  was  developed  by  Anderson 
[1973]  to  solve  for  0+  along  flux  tubes.  (0+  is  the  dominant  ion  in  the  F  region  of  the 
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ionosphere,  which  is  the  main  part  of  the  ionosphere.)  The  method  is  based  on  the  concept  of 
ambipolar  diffusion,  in  which  the  inertial  terms  of  the  momentum  equation  are  neglected  and 
charge  neutrality  (ne  =  n& )  and  parallel  flux  balance  (ne\e  -b  =  «0,v0*  b)  are  enforced.  The 
continuity  equation  for  0+  is 


The  0+  flux  may  be  obtained  from  the  electron  and  0+  momentum  equations  by  using  the  flux 
balance  condition,  neglecting  inertial  terms: 


W 


(2) 


Upon  substitution  of  Eq.  (2)  into  Eq.  (1),  one  obtains  a  linear  diffusion  equation  which  in 
LOWLAT  is  solved  using  the  Crank-Nicholson  method.  The  chief  problem  with  this  approach  is 
that  the  collision  frequency  v  „  becomes  very  small  at  high  altitudes  (where  0+  may  not  be  the 

dominant  ion)  leading  to  numerical  instabilities.  This  was  handled  in  LOWLAT  by  adding  to  the 
actual  collision  term  a  constant  so  that  the  collision  frequency  cannot  fall  below  the  constant 
value. 


When  one  extends  the  preceding  equations  to  include  additional  ions,  the  resulting 
coupled  diffusions  equations  become  nonlinear.  The  full  derivation  of  the  multi-ion  ambipolar 
diffusion  equations  is  given  in  Appendix  B.  The  impact  of  this  non-linearity  on  the  solution 
method  used  in  GTIM  had  not  been  properly  assessed  before  CPI  was  given  access  to  the  code. 
The  original  solution  method  was  simply  to  apply  the  linear  solver  to  the  three  separate  ions 
without  taking  into  account  the  additional  terms  that  appear  in  the  multi-ion  equations,  but  vanish 
when  only  one  ion  is  present.  Working  primarily  under  the  contract  F19628-96-C-0041  (for 
which  a  working,  multi-ion  version  of  GTIM  was  a  prerequisite),  we  modified  the  code  to  be 
consistent  with  the  proper  multi-ion  equations.  In  that  case,  it  does  not  appear  that  the  collision 
frequencies  need  to  be  modified  to  avoid  numerical  instabilities,  a  definite  advantage  of  this 
approach. 

As  of  August  1997,  we  were  still  studying  the  implications  of  this  situation  and  the 
behavior  of  the  code.  It  was  not  clear  that  the  solutions  being  obtained  were,  in  fact,  correct. 
Although  the  0+  densities  appeared  to  be  consistent  with  observations  and  other  models,  the  H+ 
and  He+  densities  were  behaving  properly.  We  will  report  on  the  results  of  our  investigations  and 
the  resolution  of  the  problem  in  future  reports. 
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Appendix  A 
PRISM  Updates 


Contents: 


PRISM  changes  memoranda  for  the  period  12  August  1995  through  19  August  1996. 


Memo  date 

PRISM  version 

Page 

1996  September  30 

1.6b  to  1.7 

6 

1996  November  7 

1.7  to  1.7a 

11 

1997  February  14 

1.7a  to  1.7b 

17 
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MEMORANDUM 


DATE:  30-September-1996 

TO:  Rob  Daniell 

FROM:  Lincoln  Brown 

RE:  Changes  to  PRISM  1. 6b  for  PRISM  1. 7 


The  changes  to  PRISM  1.6b  for  PRISM  1.7  focus  on  replacing  the  LLF  coefficient  set  again,  on 

improving  the  handling  of  collocated  data  in  the  real-time  adjustment,  and  on  improving  the 

merging  of  the  parameterized  model  density  profiles.  The  changes  are  summarized  as  follows: 

I.  BUG  FIX:  The  LLF  parameterized  model  coefficient  set  has  been  replaced  again. 
Because  of  a  bug  in  the  processing  of  the  LOWLAT  output,  the  magnetic  latitude  grid  in  the 
PRISM  1.6b  LLF  coefficient  files  was  defined  to  start  at  -33  °N  instead  of  the  correct  value  of 
-34°N.  This  resulted  in  an  erroneous  northward  1°  magnetic  latitude  shift  in  the  0+  densities 
from  the  LLF  parameterized  model.  The  LLF  coefficient  files  have  been  regenerated  with  the 
correct  starting  value  for  the  magnetic  latitude  grid.  Note  that  this  fix  requires  no  change  to 
the  PRISM  source  code. 

II.  BUG  FIX:  The  conversion  of  UT  from  hours  to  HHMM  in  routine  PARAM  for  the 
URSI  f0F2  model  has  been  corrected.  Previously,  the  conversion  of  UT  from  decimal  hours  to 
HHMM  resulted  in  truncation  of  the  minutes  to  zero,  e.g.  a  UT  of  1 .5  hours  was  erroneously 
converted  to  0100.  The  conversion  has  been  modified  so  that  the  minutes  are  not  truncated  to 
zero,  e.g.  a  UT  of  1 .5  hours  is  correctly  converted  to  0130. 

IE.  Based  on  work  done  by  Pat  Doherty,  we  now  recommend  that  a  27-day  running  mean  of 
Fio.7  be  used  instead  of  a  daily  value.  Pat  found  that  a  mean  F10.7  correlates  much  better  with 
the  ionosphere  than  a  daily  value,  either  because  the  solar  EUV  has  less  day-to-day  variability 
than  F10.7  or  because  the  day-to-day  variability  in  solar  EUV  is  uncorrelated  with  the  day-to- 
day  variability  in  F10.7.  The  same  applies  to  Sunspot  Number  (use  a  27-day  running  mean 
instead  of  a  daily  value),  although  we  recommend  using  F10.7  instead  of  Sunspot  Number  as 
the  solar  activity  index. 

IV.  Multiple  data  points  critically  close  to  an  output  grid  point  are  now  handled  correctly  in 
the  real-time  adjustment.  Previously,  if  one  or  more  data  points  was  within  the  critical 
distance  of  an  output  grid  point,  then  the  correction  returned  by  the  real-time  adjustment  was 
the  correction  from  the  last  data  point  within  the  critical  distance.  Now,  if  one  or  more  data 
points  are  within  the  critical  distance  of  an  output  grid  point,  then  the  correction  returned  from 
the  real-time  adjustment  is  the  average  of  corrections  from  those  data  points. 

V.  The  top-level  parameterized  model  routine  PARAM  has  undergone  substantial  changes. 
The  original  motivation  was  a  report  from  Jim  Secan  at  Northwest  Research  Associates 
illustrating  a  discontinuity  in  PRISM  1 .6b  TEC  at  the  latitude  transition  region  between  the 


mid-latitude  and  high-latitude  parameterized  models.  A  closer  examination  of  routine 
PARAM  has  yielded  the  following  improvements: 

A.  It  has  been  completely  rewritten  for  optimization  and  readability.  The  number  of 
routine  calls  has  been  reduced. 

B.  The  transition  between  the  mid-latitude  and  high-latitude  parameterized  models  is 
smoother.  Previously,  in  the  mid/high-latitude  transition  region,  the  mid-latitude  peak 
density  and  height  were  used  to  merge  the  mid-latitude  and  high-latitude  density  profiles, 
resulting  in  a  discontinuity  in  density  at  the  boundary  between  the  high-latitude  region  and 
the  mid/high-latitude  transition  region.  Now,  weighted  averages  of  mid-latitude  and  high- 
latitude  peak  densities  and  heights  are  used  in  the  merging,  resulting  in  a  smoother 
transition. 

C.  The  transition  between  the  low-latitude  and  mid-latitude  parameterized  models  is 
smoother.  Previously,  a  weighted  average  of  low-latitude  and  mid-latitude  critical 
frequencies  was  used  to  calculate  the  peak  density  for  the  merging  of  the  low-latitude  and 
mid-latitude  density  profiles.  Now,  the  peak  density  used  in  the  merging  is  calculated  from 
a  weighted  average  of  low-latitude  and  mid-latitude  peak  densities.  This  is  consistent  with 
the  profile  merging  process,  and  results  in  a  smoother  transition  between  the  low-latitude 
and  mid-latitude  parameterized  models. 

VI.  The  mid-level  parameterized  model  routines  LOW_PARAM,  MDD_PARAM,  and 
USUMODEL  have  been  rewritten  to  accommodate  changes  in  the  top-level  parameterized 
model  routine  PARAM,  to  provide  small  performance  gains  by  reducing  routine  calls  and 
eliminating  unnecessary  calculations,  and  to  improve  readability. 

VII.  The  low-level  parameterized  USU  model  routines  RECON,  FOUR_COEFF,  and 
FULL_DATA  have  been  modified  to  accommodate  changes  to  the  mid-level  parameterized 
model  routine  USUMODEL. 

VIII.  Several  routines  in  module  ENVIRON.FOR  have  been  replaced  for  compatibility  with 
PIM. 

IX.  Several  routines  have  been  modified  to  accommodate  the  new  version  of  routine 
SOLDEC. 

X.  Several  minor  changes  have  been  made  for  compatibility  with  Microsoft  FORTRAN. 
They  do  not  impact  the  results. 

The  table  below  describes  the  changes  that  I  made  to  PRISM  1 .6b  to  produce  PRISM  1 .7. 
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Module 

Program  Unit 

Description  of  Changes 

ADJ_FNH.FOR 

Subroutine  GET_FOS 

Removed  argument  ERR  from  the  call  to  INFLECTION  since  it  is  no  longer  used  by  that  routine. 
Removed  local  variable  ERR  since  it  is  no  longer  used. 

Subroutine  INFLECTION 

Removed  input  argument  ERR  since  it  is  not  used. 

Changed  declaration  of  local  variable  TEST  from  LOGICAL*!  to  LOGICAL 

Subroutine  GET_FOE 

Removed  since  it  is  no  longer  used. 

ENVIRON.FOR 

Subroutine  SOLZA 

Removed  since  it  has  been  replaced  by  routine  SOLSZA. 

Subroutine  SSOLPT 

Removed  since  it  has  been  replaced  by  routine  SOLSUB. 

Function  SOLDEC 

Replaced  by  routine  of  the  same  name. 

Function  SOLANG 

New  routine. 

Subroutine  SOLSUB 

New  routine. 

Function  SOLSZA 

New  routine. 

FMODEL.FOR 

Subroutine  FMODEL 

Removed  argument  HBOT  from  call  to  PARAM  since  it  is  not  used  by  that  routine. 

Removed  last  argument  in  call  to  PARAM  since  it  is  no  longer  used  by  that  routine. 

Removed  local  variable  HBOT  since  it  is  no  longer  used. 

GETDAT.FOR 

Subroutine  GETDAT 

Removed  argument  ESWITCH  from  call  to  MODEL_FLAGS  since  it  is  no  longer  used. 

Subroutine  MODEL_FLAGS 

Removed  assignment  of  input  argument  ESW. 

Removed  input  argument  ESW  since  it  is  no  longer  used. 

Subroutine  PREC„DATA 

Changed  declaration  of  local  variable  TEST  from  LOGICAL*!  to  LOGICAL. 

Subroutine  J4ARR 

Changed  declaration  of  local  variable  TEST  from  LOGICAL*l  to  LOGICAL. 

HL1M.FOR 

Subroutine  REGMOD 

Removed  argument  HBOT  from  call  to  PARAM  since  it  is  not  used  by  that  routine. 

Removed  last  argument  in  call  to  PARAM  since  it  is  no  longer  used  by  that  routine. 

Removed  local  variable  HBOT  since  it  is  no  longer  used. 

Replaced  call  to  SOLZA  with  call  to  SOLSZA. 

HLISM.FOR 

nmH 

Removed  local  variable  MFLAG  since  it  is  no  longer  used. 

INDIRECT.INC 

n/a 

Removed  common  block  INDIRECT  variable  ESWITCH  since  it  is  no  longer  used. 

1NIT.FOR 

Subroutine  LDITER 

Changed  declaration  of  output  argument  LAYFLG  from  LOGICAL*  1  to  LOGICAL. 

Subroutine  GETFC 

Changed  declaration  of  local  variable  ERROR  from  LOGICAL*  1  to  LOGICAL. 

Subroutine  LUDCMP 

Changed  PAUSE  statement  to  STOP  statement. 

IO_UTIL.FOR 

Subroutine  UDETID 

Added  input  argument  UT. 

Argument  DAY  is  no  longer  converted  to  a  real  number  in  call  to  SOLDEC. 

Added  argument  UT  to  call  to  SOLDEC. 

LAYFLG.INC 

n/a 

Changed  declaration  of  common  block  LAYFLG  variable  LAYFLG  from  LOGICAL*  1  to 

LOGICAL. 
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Module  I  _ Program  Unit 


LOW_PARA.FOR  Subroutine  LOW_P ARAM 


MATHJJTI.FOR  Subroutine  INTRP _ 

Subroutine  SVDCMP 


MID_PARA.FOR  Subroutine  MID_PAR  AM 


OUTPUT. FOR  I  Subroutine  FINAL 


PAR  AM.  FOR  Subroutine  PAR  AM 


PRISM. FOR 


READ_DBA.FOR 


USERJNP.FOR 


Program  PRISM 


Subroutine  READ__DBASES 


Subroutine  RDLOW 


Subroutine  LDETID 


Subroutine  READ_E 


Subroutine  EDETID 


Subroutine  READMID 


Subroutine  DETID 


Subroutine  READUSU 


Subroutine  RTA 


Subroutine  CORRECT! 


Subroutine  COR_M  AX 


Subroutine  GET  ONE„FO 


Subroutine  COR_PRO 


Subroutine  CHECK  ST AT 


Subroutine  GIVE_DATA 


Description  of  Changes  (continued) _ 


Rewrite  for  optimization  and  readability. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 
Changed  intrinsic  function  ALOG  to  its  generic  equivalent  LOG. 


Modified  logic  to  eliminate  arithmetic  IF  statements. _ _ _ ___ _ 


Changed  PAUSE  statement  to  PRINT  statement. _  — 


Rewrite  for  optimization  and  readability. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 
Changed  intrinsic  function  ALOG  to  its  generic  equivalent  LOG. 


Changed  declaration  of  local  variable  TEST  from  LOGICAL*!  to  LOGICAL.  — 


Rewrite  for  optimization  and  readability. 

Fixed  a  bug  in  the  conversion  of  UT  from  hours  to  HHMM  for  the  call  to  F2URS1. 

In  the  mid/high-latitude  transition  region,  weighted  averages  of  mid-latitude  and  high-latitude  peak 
densities  and  heights  are  now  used  in  the  merging  of  profiles  instead  of  the  mid-latitude  peak  density 
and  height. 

In  the  low/mid-latitude  transition  region,  the  peak  density  used  to  merge  profiles  is  calculated  from  a 
weighted  average  of  low-latitude  and  mid-latitude  peak  density  instead  of  a  weighted  average  of  low- 
latitude  and  mid-latitude  critical  frequency,  to  be  consistent  with  the  profile  merging  process. 
Simplified  logic  since  LAYR  flag  is  no  longer  used. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 

Removed  input  argument  HBOT  since  it  is  not  used. 

Removed  last  argument  from  calls  to  USUMODEL  since  it  is  no  longer  used. 

A  single  call  to  USUMODEL  is  now  used  instead  of  two  calls. 

Removed  last  argument  from  calls  to  MID__PARAM  since  it  is  no  longer  used. 

A  single  call  to  MID_PARAM  is  now  used  instead  of  two  calls. 

Removed  last  argument  from  calls  to  LOW_PARAM  since  it  is  no  longer  used. 

A  single  call  to  LOWJPARAM  is  now  used  instead  of  two  calls. 

Removed  call  to  GET_FOE  since  it  is  no  longer  used. 

Removed  logic  involving  common  block  INDIRECT  variable  ESWITCH  since  it  is  no  longer  used. 
Changed  intrinsic  function  ALOG  to  its  generic  equivalent  LOG. 


Updated  the  version  number  and  version  date. _ 


Added  argument  UT  to  call  to  RDLOW. _ _ _ _ _ 


Added  input  argument  UT. 

Added  argument  UT  to  call  to  LDETID.  _ _ _ 


Added  input  argument  UT. 

Argument  DAY  is  no  longer  converted  to  a  real  number  in  call  to  SOLDEC. 

Added  argument  UT  to  call  to  SOLDEC.  _ _ _ 


Added  argument  UT  to  call  to  EDETID.  _ _ _ _ 


Added  input  argument  UT. 

Argument  DAY  is  no  longer  converted  to  a  real  number  in  call  to  SOLDEC. 

Added  argument  UT  to  call  to  SOLDEC. _ _ _ 


Added  argument  UT  to  call  to  DETID.  _ _ _ 


Added  input  argument  UT. 

Argument  DAY  is  no  longer  converted  to  a  real  number  in  call  to  SOLDEC. 

Added  argument  UT  to  call  to  SOLDEC.  _  _ 


Added  argument  UT  to  calls  to  UDETID.  _ 


Changed  declaration  of  local  variable  IERR  from  LOGICAL*!  to  LOGICAL. _ 

Rewrite  for  optimization  and  readability. 

If  one  or  more  data  points  is  within  the  critical  distance  of  the  point  of  interest,  then  the  correction  at 
the  point  of  interest  is  now  calculated  as  the  average  of  corrections  from  those  data  points  instead  of 
the  correction  of  the  last  data  point  within  the  critical  distance  being  used.  _ 


Changed  declaration  of  output  argument  IERR  from  LOGICAL*!  to  LOGICAL. _ 


Changed  declaration  of  output  argument  IERR  from  LOGICAL*!  to  LOGICAL. _ 


Changed  declaration  of  output  argument  IERR  from  LOGICAL*!  to  LOGICAL. 


Changed  declaration  of  input  argument  TEST  from  LOGICAL*!  to  LOGICAL.  _ 


Changed  declaration  of  local  variable  TEST  from  LOGICAL*  1  to  LOGICAL. _ _ 


Module 

HNBMI JjMJ  Will 

Description  of  Changes  (< continued) 

USUMODEL.FOR 

Subroutine  USUMODEL 

Rewrite  for  optimization  and  readability. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 

Changed  intrinsic  function  ALOG  to  its  generic  equivalent  LOG. 

Subroutine  RECON 

Replaced  input  argument  XNMLAT  with  input  arguments  SMLATE  and  SMLATF. 

Replaced  argument  XNMLAT  with  argument  SMLATE  and  SMLATF  in  call  to  FOUR_COEFF. 
Removed  argument  LAYR  from  calls  to  FOUR_COEFF  and  FULL_DATA. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 

Subroutine  FOUR_COEFF 

Replaced  input  argument  XNMLAT  with  input  arguments  SMLATE  and  SMLATF. 

Fourier  coefficients  are  now  calculated  for  all  three  ions. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 

Subroutine  FULL_DATA 

Altitude  profiles  are  now  calculated  for  all  three  ions. 

Removed  input  argument  LAYR  since  it  is  no  longer  used. 

The  figure  below  illustrates  the  improvement  in  the  merging  of  the  regional  parameterized 
models.  It  shows  a  profile  of  f0F2  vs.  geographic  latitude.  Note  the  discontinuity  in  PRISM  1 .6b 
foF2  in  the  mid/high-latitude  transition  region  and  the  smoother  transition  in  the  PRISM  1.7 
profile.  The  differences  at  low  latitudes  are  due  to  the  change  in  the  low/mid-latitude  transition 
method  and  the  corrected  LLF  parameterized  model  coefficients. 


test. out 

Year:  1996 
Day:  80 
UT:  12.00  hr 

F10.7:  193.0 
SSN:  150.0 
Kp:  3.0 

Lon:  150.0  °E  geographic 
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MEMORANDUM 


DATE:  7 -November- 1996 

TO:  Rob  Daniell 

FROM:  Lincoln  Brown 

RE:  Changes  to  PRISM  1. 7 for  PRISMUa _ ^ 

The  changes  to  PRISM  1.7  for  PRISM  1.7a  focus  on  improvements  to  the  real-time  data 

ingestion.  The  changes  are  summarized  as  follows: 

I  BUG  FIX:  The  presence  of  all  five  possible  DISS  real-time  data  files  no  longer 
causes  an  array-out-of-bounds  error  in  the  DISS  data  ingestion.  Previously,  if  root  names 
for  all  five  DISS  data  files  were  given  in  the  PATH_NAM.TXT  file,  then  the  DISS  data 
ingestion  would  attempt  to  access  a  nonexistent  sixth  element  of  the  array  that  holds  the 
DISS  data  file  root  names.  This  problem  has  been  eliminated  due  to  the  rewrite  of  the 
DISS  data  ingestion  (see  item  5). 

II.  BUG  FIX:  The  presence  of  all  five  possible  IMS  real-time  data  files  no  longer  causes 
an  array-out-of-bounds  error  in  the  IMS  data  ingestion.  Previously,  if  root  names  for  all 
five  IMS  data  files  were  given  in  the  PATH_NAM.TXT  file,  then  the  IMS  data  ingestion 
would  attempt  to  access  a  nonexistent  sixth  element  of  the  array  that  holds  the  IMS  data 
file  root  names.  This  problem  has  been  eliminated  due  to  the  rewrite  of  the  IMS  data 
ingestion  (see  item  6). 

III.  BUG  FIX:  The  presence  of  all  eight  possible  DMSP  real-time  data  files  no  longer 
causes  an  array-out-of-bounds  error  in  the  DMSP  data  ingestion.  Previously,  if  root  names 
for  all  eight  DMSP  data  files  were  given  in  the  PATH_NAM.TXT  file,  then  the  DMSP 
data  ingestion  would  attempt  to  access  a  nonexistent  sixth  element  of  the  array  that  holds 
the  DMSP  data  file  root  names.  This  problem  has  been  eliminated  due  to  the  rewrite  of 
the  DMSP  data  ingestion  (see  item  7). 

IV.  BUG  FIX:  Large  amounts  of  DISS  foF2  real-time  data  and/or  large  amounts  of  IMS 
real-time  data  no  longer  cause  an  array-out-of-bounds  error  in  the  midlatitude  correction 
algorithm.  Previously,  array  allocation  in  the  midlatitude  correction  algorithm  for 
corrections  for  these  kinds  of  data  was  inconsistent  with  limits  defined  in  the  real-time 
data  ingestion.  This  problem  was  reported  by  David  Coxwell  at  Air  Force  Institute  of 
Technology.  This  problem  has  been  eliminated  due  to  the  rewrite  of  the  DISS  and  IMS 
data  ingestion  (see  items  5  and  6)  and  the  rewrite  of  the  midlatitude  correction  algorithm 
(see  item  9). 

V.  The  DISS  real-time  data  ingestion  has  been  rewritten  for  speed  and  readability.  A 
single  pass  through  the  DISS  data  is  now  required  instead  of  two  passes.  Up  to  850  DISS 
data  records  can  now  be  accepted  (i.e.  meet  criteria  for  use  in  the  real-time  adjustment) 
from  all  DISS  data  files  combined,  a  reduction  from  the  previous  limit  of  7000  accepted 


records  per  DISS  data  file.  The  limit  of  850  accepted  DISS  data  records  is  based  on  50 
DISS  stations  and  a  15  minute  time  resolution  in  a  [-2,2]  hour  time  window  relative  to  the 
time  of  the  run. 

VI.  The  IMS  data  ingestion  has  been  rewritten  for  speed  and  readability.  A  single  pass 
through  the  IMS  data  is  now  required  instead  of  two  passes.  Up  to  21600  IMS  data 
records  can  now  be  accepted  (i.e.  meet  criteria  for  use  in  the  real-time  adjustment)  from  all 
IMS  data  files  combined,  a  reduction  from  the  previous  limit  of  7000  accepted  records  per 
IMS  data  file.  The  limit  of  21600  accepted  IMS  data  records  is  based  on  200  IMS  stations 
reporting  12  simultaneous  measurements  and  a  time  resolution  of  15  minutes  in  a  [-2,2] 
hour  time  window  relative  to  the  time  of  the  run. 

VII.  The  DMSP  data  ingestion  has  been  rewritten  for  speed  and  readability.  A  single  pass 
through  the  DMSP  data  is  now  required  instead  of  three  passes.  Up  to  24  sets  of  each  of 
the  three  kinds  of  DMSP  data  are  now  accepted  by  PRISM  from  all  DMSP  data  files 
combined,  a  reduction  from  the  previous  limit  of  24  sets  per  DMSP  data  file.  The  limit  of 
24  sets  is  based  on  8  satellites,  an  orbital  period  of  6060  seconds,  and  a  [-2,2]  hour  time 
window  relative  to  the  time  of  the  run.  Up  to  1451  SSIES  ION  DRIFT  data  records  can  be 
accepted  (i.e.  meet  criteria  for  use  in  determining  equatorward  trough  boundaries)  per  set, 
a  limit  which  has  not  changed.  The  limit  of  1451  accepted  SSIES  ION  DRIFT  data 
records  is  based  on  a  5  second  time  resolution  in  a  [-2,0]  hour  time  window  relative  to  the 
time  of  the  run,  with  10  extra  records  allowed  for  orbital  overlap.  Up  to  3848  SSIES  IN 
SITU  PLASMA  data  records  can  now  be  accepted  (i.e.  meet  criteria  for  use  in  the  real¬ 
time  adjustment)  from  all  DMSP  data  files  combined,  a  reduction  from  the  previous  limit 
of  1451  accepted  records  per  DMSP  data  file.  The  limit  of  3848  accepted  SSIES  IN  SITU 
PLASMA  data  records  is  based  on  8  satellites  and  a  time  resolution  of  30  seconds  in  a  [- 
2,2]  hour  time  window  relative  to  the  time  of  the  ran.  Up  to  72 1 1  SSJ/4  data  records  can 
be  accepted  (i.e.  meet  criteria  for  use  in  determining  auroral  oval  boundaries  and  for  use  in 
the  real-time  adjustment)  per  set,  a  limit  which  has  not  changed.  The  limit  of  7211 
accepted  SSJ/4  data  records  is  based  on  a  time  resolution  of  1  second  in  a  [-2,0]  hour  time 
windows  relative  to  the  time  of  the  run,  with  1 0  extra  records  allowed  for  orbital  overlap. 

VIE.  The  output  station  list  determination  has  been  rewritten  for  speed  and  readability.  A 
single  pass  through  the  output  station  list  is  now  required  instead  of  two  passes.  The 
number  of  output  stations  is  now  unlimited,  a  change  from  the  previous  limit  of  7000 
stations. 

IX.  The  midlatitude  correction  algorithm  has  been  rewritten  for  the  following  reasons: 

A.  To  optimize  for  speed. 

B.  To  prevent  extraneous  zero  corrections  from  being  included  in  the  real-time 
adjustment.  The  extraneous  corrections  were  not  based  on  valid  data,  and,  although 
zero,  influenced  the  correction  fields  and  increased  run  time. 

C.  To  remove  coding  for  BOTTOMSIDE  data  since  that  data  type  is  no  longer  considered  by 
PRISM. 

D.  To  remove  redundant  and  commented-out  coding. 

E.  To  improve  readability. 

X.  Limits  on  the  ingestion  of  real-time  data  are  now  globally  defined  in  new  INCLUDE 
file  rtdlimit.inc. 
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XI.  The  HLE  model  has  been  removed  since  it  is  no  longer  used.  All  components  of 
PRISM  upon  which  HLE  solely  depends  have  also  been  removed,  such  as  the  MSIS-86 
neutral  atmosphere  model  and  a  number  of  INCLUDE  files.  The  HLE  data  files 
CHEM.FIL,  TIMING.FIL,  and  SOLAR.FIL  are  no  longer  used.  The  path  to  the  HLE  files 
is  still  read  from  the  PATH_NAM.TXT  file  for  compatibility,  but  it  is  not  used. 

XD.  Coding  for  the  BOTTOMSIDE  data  type  has  been  removed  since  it  is  no  longer 
considered  in  PRISM. 

XIII.  INCLUDE  file  dlta.inc  has  been  removed  since  it  is  not  used. 

XIV.  A  number  of  unused  PARAMETERS  and  variables  have  been  removed. 

The  performance  gain  due  to  these  changes  is  hard  to  predict  since  it  depends  on  the  amount  and 
character  of  the  real-time  data.  For  many  of  the  standard  test  cases  provided  with  delivery,  run¬ 
time  has  been  decreased  by  about  15%.  The  memory  requirement  has  been  reduced  by  about  0.5 
MB,  a  7%  decrease. 

The  table  below  describes  the  changes  that  I  made  to  PRISM  1 .7  to  produce  PRISM  1 .7a. 
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Module 


ALPHAN.1NC 


CHEM.INC 


DLTA.INC 


ESPECT.INC 


n/a 


n/a 


n/a 


n/a 


Subroutine  BNDRY 


Subroutine  NEUATM 


Subroutine  GTS5 


Function  DENSS 


Function  GLOBES 


Function  GLOB5L 


Function  DNET 


Function  CCOR 


Block  Data  PRMDTD 


Function  SOLANG 


Subroutine  SOLSUB 


Function  SOLSZA 


Subroutine  GETDAT 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


DO  DIR 


DOJONO 


DO  TEC 


DO_DMSP 


INRTD 


1NDISS 


INIMS 


INDMSP 


DETOSL 


PREC_DATA 


Subroutine  LOADDIR 


HLEMODEL.FOR  n/a 


HLIM.FOR  Subroutine  REGMOD 


Description  of  Changes 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Added  INCLUDE  statement  for  INCLUDE  File  rtdlimit.inc. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Replaced  call  to  routine  DOJDIR  with  calls  to  routine  INRTD  and  routine  DETOSL. 
Removed  output  parameter  SUSI  since  it  is  not  used. 

Removed  PARAMETERS  NSTA  and  MAXTYPE  since  they  are  no  longer  used. 
Removed  local  variable  NDAT  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used  (it  has  been  replaced  by  routine  INRTD). 


Removed  since  it  is  no  longer  used  (it  has  been  replaced  by  routines  INDISS  and  DETOSL). 


Removed  since  it  is  no  longer  used  (it  has  been  replaced  by  routine  INIMS). 


Removed  since  it  is  no  longer  used  (it  has  been  replaced  by  routine  INDMSP). 


New  routine. 


New  routine. 


New  routine. 


New  routine. 


New  routine. 


Output  parameter  NR  EC  is  now  initialized  to  zero. 

PRECIPITATION  records  are  now  written  to  data  file  at  the  current  file  position  instead  of  the  end  of 
the  file. 

A  PRECIPITATION  data  header  is  no  longer  written  to  the  data  file. 

Standard  deviation  placeholders  are  no  longer  included  in  the  data  record  written  to  the  data  file. 
INCLUDE  statements  for  INCLUDE  files  dpath.inc  and  directed. inc  have  been  removed  since  they 
are  no  longer  needed. 

Removed  local  variable  DUMMY  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 


Removed  coding  for  HLE  model  since  it  is  no  longer  used. 

Removed  input  parameters  F10ANA,  RF10NA,  and  APNA  since  they  are  no  longer  used. 
Removed  local  variables  SZA,  HMED,  and  ZNEU  since  they  are  no  longer  used. 

Removed  PARAMETER  NZNEU  since  it  is  no  longer  used. 

Removed  coding  for  BOTTOMSIDE  data  since  the  BOTTOMSIDE  data  type  is  no  longer  used. 


HLISM.FOR 

Subroutine  HLISM 

Removed  arguments  RF10P7,  F10P7A,  and  AP  from  call  to  routine  MATRIX  since  they  are  no 
longer  used  by  that  routine. 

Removed  input  parameters  RF1 0P7,  F10P7A,  and  AP  since  they  are  no  longer  used. 

Subroutine  MATRIX 

Removed  arguments  F10P7A,  RF10P7,  and  AP  from  calls  to  routine  REGMOD  since  they  are  no 
longer  used  by  that  routine. 

Removed  input  parameters  RF10P7,  F10P7A,  and  AP  since  they  are  no  longer  used. 

INIT.FOR 

Subroutine  INITPR 

Added  INCLUDE  statement  for  INCLUDE  file  rtdlimit.inc. 

PARAMETER  MAXORB  has  been  renamed  to  MORBIT. 

INT.INC 

n/a 

Removed  since  it  is  no  longer  used. 

LOG  1C  UN  I.  INC 

n/a 

Removed  PARAMETERS  LUIDL,  LUNATM,  LUNUMB,  LURATE,  LURI,  LURO.  LUSGP,  and 
LUTEMP  since  they  are  not  used. 

Removed  PARAMETERS  LUCHEM,  LUSOLR,  and  LUT1ME  since  they  are  no  longer  used. 

LT.INC 

n/a 

Removed  since  it  is  no  longer  used. 

Module 


MATHJJTI.FOR  Subroutine  RTB1S 
Subroutine  1NTRP 
Subroutine  CHMRRC 


MATH4_CO.INC _ _ n/a_ 


M1DLAT.FOR  Subroutine  MIDLAT 


Subroutine  PER_ARR 


NATMOS.INC 


NEWFIT.INC 


Block  Data  PREC 


Program  PRISM 


PSPECT.INC 


READ_DBA.FOR  Subroutine  DBASES 

Subroutine  READ__DBASES 

Subroutine  RDF1LES _ 

Subroutine  RPC  HEM _ 

Subroutine  RDT1ME _ 

Subroutine  RDSOLR  _ 


Description  of  Changes  ( continued) 


Removed  since  it  is  no  longer  used. _ _ _ 


Removed  since  it  is  no  longer  used. _ _ _ _ — . 


Removed  since  it  is  no  longer  used.  _ 


OUTPUT.FOR  I  Subroutine  GRID_OUTPUT 


Subroutine  STATION_OUTPUT 


PHANTOM. FOR  I  Subroutine  PH ANTM 


PHYS4  CO.INC 


PRECIP.FOR 


PRECIP.INC 


PRISM. FOR 


Rewrite  for  the  following  reasons: 

1 .  To  optimize  calculation  of  corrections. 

2.  To  prevent  extraneous  zero  corrections  from  being  included  in  the  real-time  adjustment.  The 
extraneous  corrections  were  not  based  on  valid  data,  and  although  zero,  did  influence  the 
correction  fields,  and  increased  run  time. 

3.  To  remove  coding  for  BOTTOMS1DE  data  since  that  data  type  is  no  longer  used. 

4.  To  remove  redundant  and  commented-out  coding. 

5.  To  improve  readability.  .  , 

Removed  arguments  F10P7A,  RFI0P7,  and  AP  from  calls  to  routine  REGMOD  since  they  are  no 

longer  used  by  that  routine. _ _ _ _ _ 

Removed  since  it  is  no  longer  used 


Removed  since  it  is  no  longer  used. 


Removed  since  it  is  no  longer  used. 

Combined  common  blocks  INWFT  and  NEWFT  into  common  block  NEWFIT 
Removed  variables  NHFM,  NFFM,  NHEM,  NFEM,  NB1 M,  NBI D,  NB2M,  NB2D,NHTM,  NNTM, 
FFDATA,  EFDATA,  FHDATA,  EHDATA,  BIDATA,  B2DATA,  NTDATA,  HTDATA,  BT1LAT, 
BTILOn!  BT2LAT,  BT2LON,  BIC,  and  B2C  from  common  block  NEWFIT  since  they  are  no 

Replaced  PARAMETERS  NMAX  and  NNMAX  with  PARAMETERS  MFFD,  MHFD,  MFED, 
MHED,  MNTD,  and  MHTD.  The  new  PARAMETERS  are  dependent  on  PARAMETERS  in 

INCLUDE  file  rtdlimit.inc.  - 

Removed  arguments  F10P7A,  RF10P7,  and  AP  from  call  to  routine  REGMOD  since  they  are  no 
longer  used  by  that  routine. 

Removed  INCLUDE  statement  for  INCLUDE  file  dlta.inc  since  it  is  not  used. _ _ _ 

Removed  arguments  F10P7A,  RFI0P7,  and  AP  from  call  to  routine  REGMOD  since  they  are  no 
longer  used^Jha^outine^^^^^||([^^[|||^|^^^^^^^^^ll<^^^l^^^^^^^^^,^^^w. 

Removed  arguments  FI0P7A,  RFI0P7,  and  AP  from  calls  to  routine  REGMOD  since  they  are  no 
longer  used  by  that  routine. 

Removed  INCLUDE  statement  for  INCLUDE  file  indirect.inc  since  it  is  no  longer  needed, _ 


Removed  since  it  is  no  longer  used. _  ... 


Added  INCLUDE  statement  for  INCLUDE  file  rtdlimit.inc. 


Added  INCLUDE  statement  for  INCLUDE  file  rtdlimit.inc. 

PARAMETER  MAXORB  has  been  renamed  to  MORB1T. _ _ _ 

Removed  argument  SUSI  from  call  to  routine  GETDAT  since  it  is  not  used. 

Removed  arguments  RF10P7,  F10P7A,  and  AP  from  calls  to  routine  HLISM  since  they  are  no  longer 
used  by  that  routine. 

Removed  INCLUDE  statement  for  INCLUDE  file  dlta.inc  since  it  is  not  used. 

Updated  the  version  number  and  version  date. 


Removed  since  it  is  no  longer  used. 

Removed  argument  PLHE(J9:J10)  from  call  to  routine  READ_DBASES  since  it  is  no  longer  used  by 
that  routine. 

Removed  local  variables  J9  and  JO  since  they  are  no  longer  used. _ _ _ _ 


Removed  call  to  routine  RDFILES  since  it  is  no  longer  needed. 

Removed  input  parameter  HLEPATH  since  it  is  no  longer  used.  — 


Removed  since  it  is  no  longer  used. _ _ 


Removed  since  it  is  no  longer  used. _ _ _ _ _ — - 


Removed  since  it  is  no  longer  used. _ _ _ _ _ 


Removed  since  it  is  no  longer  used. 


Module 


Module 

HBniRIHil 

Description  of  Changes  (< continued) 

RTA.FOR 

Subroutine  RTA 

Added  INCLUDE  statement  for  INCLUDE  file  rtdlimit.inc. 

Removed  arguments  FFDATA,  EFDATA,  FHDATA,  EHDATA,  NTDATA,  and  HTDATA  from 
calls  to  routine  CORRECT  1  since  they  are  no  longer  used  by  that  routine. 

Removed  arguments  NMAX  and  NNMAX  from  calls  to  routine  CORRECT  1  since  they  are  no 
longer  used  by  that  routine. 

Removed  coding  for  BOTTOMSIDE  data  since  the  BOTTOMSIDE  data  type  is  no  longer  used. 
Removed  INCLUDE  statement  for  INCLUDE  file  dlta.inc  since  it  is  not  used. 

Subroutine  CORRECT  I 

Validity  checking  is  no  longer  needed  since  all  data  points  reaching  this  routine  are  to  be  used. 
Removed  input  argument  MDAT  since  it  is  no  longer  used. 

Removed  input  argument  LDAT  since  it  is  no  longer  used. 

RTDLIMIT.INC 

n/a 

New  file. 

SOLAR.  INC 

n/a 

Removed  since  it  is  no  longer  used. 

SPECIE. INC 

n/a 

Removed  since  it  is  no  longer  used. 

STRINGS.FOR 

Subroutine  STRTRM 

New  routine. 

TIME. INC 

n/a 

Removed  since  it  is  no  longer  used. 

MEMORANDUM 


DATE:  1 4-February- 1997 

TO:  Rob  Daniell 

FROM:  Lincoln  Brown 

RE:  Changes  to  PRISM  1. 7 a  for  PRISM  L7b _  ^ 

The  changes  to  PRISM  1.7a  for  PRISM  1.7b  focus  on  miscellaneous  bug  fixes.  The  changes  are 
summarized  as  follows: 

1  BUG  FIX:  The  maximum  allowed  number  of  real-time  IMS  data  points  that  PRISM  can 
accept  has  been  changed  from  21,600  to  40,800  in  order  to  be  consistent  with  the  I/O 
Specification  (“200  IMS  stations  reporting  12  simultaneous  measurements  and  a  time 
resolution  of  15  minutes  in  a  +/-  2  hour  time  window  centered  on  the  date  and  UT  of  the 
run”,  or  200- 12-[1  +4-60/1 5]=40, 800).  The  previous  value  of  21,600  was  consistent  with  a 
time  resolution  of  30  minutes  (200  12  [l+4-60/30]=21,600). 

2.  BUG  FIX:  Array  allocation  for  real-time  data  has  been  increased  to  allow  for  phantom 
data.  Previously,  the  presence  of  the  maximum  allowed  number  of  real-time  data  points  in 
combination  with  phantom  data  points  could  cause  an  array-out-of-bounds  error  in  the 
calculation  of  the  low/midlatitude  correction  fields. 

3  BUG  FIX:  In  the  determination  of  the  output  station  list  from  accepted  real-time  DISS 
data,  the  search  of  the  direct  data  file  for  the  IONOSONDE  data  section  has  been  modified 
to  prevent  a  read-past-end-of-file  error.  In  subroutine  DETOSL  in  module 
GETDAT.FOR,  variable  NSTOUT  is  now  reinitialized  to  zero  after  a  non-IONOSONDE 
data  section  in  the  direct  data  file  is  skipped  over  during  the  search.  Previously,  if  the 
direct  data  file  contained  one  or  more  non-IONOSONDE  data  sections,  but  contained  no 
IONOSONDE  data  section,  the  value  of  NSTOUT  from  the  last  non-IONOSONDE  data 
section  would  be  used  to  attempt  to  read  an  IONOSONDE  data  section  past  the  end  of  the 
direct  data  file,  resulting  in  a  run-time  error. 

4.  BUG  FIX:  In  the  determination  of  the  output  station  list  from  accepted  real-time  DISS 
data,  the  corrected  geomagnetic  local  times  of  the  output  stations  now  correspond  to  the 
day  of  the  year  and  Universal  Time  of  the  run.  Previously,  the  corrected  geomagnetic 
local  times  of  the  DISS  data  records  were  used  for  the  output  stations,  resulting  in  as  much 
as  a  2  hour  error  in  corrected  geomagnetic  local  time. 

PRISM’s  memory  requirement  has  been  increased  by  about  0.25  MB,  a  3%  increase.  Its  run-time 
will  be  unaffected  unless  the  output  station  list  is  taken  from  a  very  large  set  of  accepted  DISS 
data. 

The  table  below  describes  the  changes  that  I  made  to  PRISM  1 .7a  to  produce  PRISM  1 .7b. 
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Appendix  B 

Derivation  of  the  Multi-ion  Ambipolar  Diffusion  Equation  for  GTIM 
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20  November  1996 
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(Internal  CPI  Document) 

The  material  in  this  Appendix  was  largely  developed  under 
SBIR  contract  F19628-96-C0041  and  appeared  as  Appendix  B  of 
PL-TR-97-2118  (“A  New,  Improved  Ionospheric  Correction 
Algorithm  for  Single  Frequency  GPS  Receivers”). 
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Derivation  of  the  Multi-ion  Ambipolar  Diffusion  Equation  for  GTIM 


Robert  E.  Daniell,  Jr.  and  Robert  W.  Simon 
Computational  Physics,  Inc. 

20  November  1996 
Revised  29  August  1997 


0.  Notation 

Although  we  have  attempted  to  stay  close  to  the  notation  of  Anderson  [1973],  we  found  it 
necessary  to  change  the  notation  in  some  cases.  Here  is  a  summary. 

Subscripts: 

a  denotes  charged  species,  ions  and  electrons 
s  denotes  any  species,  charged  or  neutral 

ij  denote  positive  ions 

m,  n  denote  neutral  species 
e  denotes  electrons 

Scalar  quantities 

n  denotes  concentration  (number  density)  (m'3) 

q  denotes  charge  (C) 

m  denotes  mass  (kg) 

n  denotes  volume  production  rate  (m"  s' )  due  to  photoionization  or  chemical  reaction 
A  denotes  volume  loss  frequency  due  to  chemistry  (s'1) 
v  denotes  collision  frequency  (s'1) 

03  -  qB/m  denotes  the  gyrofrequency  (rad  s'1) 

T  denotes  temperature  (K) 

denotes  the  Euler  potentials 

v  =  v  •  b  denotes  the  velocity  component  parallel  to  the  magnetic  field  (corresponds  to 
Dave’s  and  f/j 

k  denotes  Boltzmann’s  constant  (1.3807  x  10“23  J  K'1) 
e  denotes  the  magnitude  of  the  charge  on  an  electron  (1.6022  x  10~19  C) 

Vector  quantities 

v  denotes  velocity  (ms'1) 

E  denotes  electric  field  (V  m'1) 

B  denotes  magnetic  induction  (T) 

b  =  B/B  denotes  a  unit  vector  parallel  to  B  (similar  to  Dave’s  it ,  but  opposite  in  sign.) 
g  denotes  gravitational  acceleration  (m  s'2) 

Tensor  quantities 

P  denotes  the  pressure  tensor  (stress  tensor)  (Pa  =  N  m'2) 
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1.  Basic  Equations 


The  basic  equation  used  in  ionospheric  modeling  is  the  continuity  equation: 

^-+V-(nava)  =  na-Xentt  (U> 

ct 

The  terms  on  the  right  hand  side  of  the  equation  represent  the  volume  photochemical  production 
{na)  and  loss  (Aa)  rates.  The  subscript  a  denotes  any  charged  species,  i.e.,  both  electrons  and 
ions.  (In  this  and  all  subsequent  equations,  bold  serif  denotes  vectors  and  matrices,  while  bold 
sans  serif  denotes  tensors.)  In  order  to  solve  the  continuity  equation  for  na(\,t),  we  need  an 
expression  for  va(x,/),  which  we  may  obtain  from  the  momentum  conservation  equation: 


— (namaxa)  +  V •{nama\aya)  =  namag  +  naqa(E  +  \a  xB)-V-Pa  ^namava^ a  v.s)  O-2) 

dt  * 


The  last  term  on  the  right  hand  side  is  the  rate  of  change  of  momentum  carried  by  species  a  as  a 
result  of  collisions  with  other  particles,  denoted  by  subscript  s.  (Note  that  a  ranges  over  charged 
species  i  and  e,  while  5  ranges  over  both  charged  species  and  neutral  species  n.  At  this  point,  if 
considering  all  momentum  sources  and  sinks  (including  photochemical  processes),  the  proper 
step  would  be  to  use  the  continuity  equation  to  transform  the  left-hand  side  into  the  substantial 
derivative  of  v„.  However,  for  ionospheric  modeling  purposes,  these  additional  terms  are 
negligible  compared  to  the  effects  of  collisions,  so  at  this  point  we  simply  neglect  inertial  terms. 
That  is,  we  assume  that  the  ionospheric  plasma  is  in  approximate  force  balance. 

«a^g  +  «a^(E  +  VaxB)-V-Pff-X«^^K-V,)  =  0 


which  may  be  rearranged  to  produce 


vana\a  - (oa (na\a  x  b)  =  nag  +  ^ - E-—S7-P a+na^ ^ 


where  v  =  Yv  ,  <y  =  (the  gyrofrequency),  and  b  =  B  /  5  is  a  unit  vector  in  the  direction 

°  7fa  ^  °  ma 

of  the  magnetic  field.  For  GTIM  we  require  only  the  component  parallel  to  the  magnetic  field. 


'«»«*>«  =  na% *' b  +  ■  b - -j-' V •  Pa  •  b  +  na Yj Vk1 v* 


If  we  now  assume  that  the  pressure  tensor  P„  is  isotropic  and  represented  by  the  ideal  gas 
law,  then  V-Pa  =  V(nakTa).  From  this  we  may  derive  the  usual  expression  for  the  flux  in  terms 
of  a  diffusion  coefficient  and  the  density  gradient. 
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(1.6) 


«a^=- 


~V(nJa)-b  +  nag.b  +  ^E-b  +  naZvapu 

Wl„  Win.  R 


Note  that  because  we  have  stopped  with  the  momentum  equation  and  neglected  the  energy  and 
heat  flow  equations,  we  do  not  obtain  the  “thermal  diffusion”  terms  of  St. -Maurice  and  Schunk 
[1976],  which  also  appear  in  Schunk  and  Nagy  [1980]  and  Young  et  al.  [1980a,  1980b]. 


2.  General  Equations  for  Multi-species  Ionospheric  Diffusion  and  Chemistry 

If  the  scale  length  of  the  plasma  is  much  greater  than  a  Debye  Length,  then  charge 

imbalances  will  not  persist  much  longer  than  r  =  — .  This  means  that  as  long  as  we  do  not 

p  0) 

p 

wish  to  model  plasma  waves,  we  may  impose  both  the  charge  neutrality  condition 


and  the  flux  balance  condition 


ne  = 


ne»e  =  2>, 


(2.1) 


(2.2) 


where  e  denotes  electrons  and  i  denotes  ions.  The  flux  balance  condition  is  the  necessary 
condition  for  ambipolar  diffusion,  and  allows  us  to  eliminate  E  from  the  force  balance  equations 
by  writing  the  electron  equation  as 


Eb  = 

e 


^ve{neue)  +  -V(neTe)-b-meg-b-me'£ve/,u/j 

tl0 


This  result  may  be  substituted  into  the  equation  for  ion  species  i  to  yield 
if  k 

- b-V(«,.7;)  +  w,  g  •  b  +  w,.  £  visvs 


mi 


mfi 


wi„ 


Yjnj°j  + — •> •  V (neTe )-meg-b-mey£  vesvs 

m  n 


V  J  J 


(2.3) 


where  j,  like  i,  denotes  only  positive  ions.  Rearranging  and  collecting  terms  results  in 
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w=i{-^[b'v(^',+f:^b'v("'r')]+++^}b) 


+  — 1«/Z  Vis  + 

Vi  si 


melLv 

_rn^qLnL 

(  v 

1/  7  Yl  l)  . 

yeLjnJUJ 

mt  e 

m i  e  ne 

\  J  A 

where  the  sum  over  j  is  over  all  ion  species.  To  simplify  the  summation,  we  assume 
v  =  v  =o.  The  above  expression  may  be  simplified  somewhat  if  we  assume  that  all  ion 

species  i  (or  j)  are  singly  charged  positive  ions  so  that  qt  -  e : 

n,v,  =  -{- —  [b ■  V(n,T< )  +  b ■  V(n,T,) 1  +  n Jl  +  ^lg  b| 

V,  ni:  ne  .  mi .  J 

1  (2.5) 

1  tl;  V1 

+  -  «/Z  -  vZW 

[  5  L  mi  j  m‘  n<!  j 

j0  produce  a  useful  expression,  we  must  collect  all  terms  involving  the  flux  of  the  zth 
species  on  the  left  hand  side  of  the  equation.  To  do  so,  however,  we  must  first  examine  the 


collision  terms: 


/  \  (  m  ^ 

«,Y  =  ni  —  VeiVi  +  niVieVe  +  ”,  Z 

mi  J  nil  s*eA  mi  J 


The  first  term  may  be  rewritten 


m„  w, 

—  Vei°i  =  ~  Vieni°i 
m.  n„ 


where  we  have  used  the  relation  njtnjV  ie  =  nemev  ei  [Schunk  and  Nagy,  1980].  Similarly,  the 
second  term  may  be  rewritten 


niVie°e 


”iVieVe 


ft  •  fl:  V1 

—  —L.  y  n  v  —  —  v  7  n  v 

vieneue  M  ie  La  J  J 

ne  ne  j 

n  n  ■r-i 

=  ~  VleniVi  +  Vie Z  n.i°i 


V  J*i 


Thus,  Eq.  (2.6)  may  be  rewritten  as 


,  n  \ 

«/Z  v«  +  —  v*  k  =  2—  KAV,-  +  7-  vteZwy°y  +  2 

s  ^  Wlj  J  ne  j*i  s*ej\  mi  J 
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This  is  to  be  combined  with  the  other  collision  term  so  that 


+- 


ma 


- 


m; 


a%!h 

mt  ne 


f 

V  J 

+r-^Sw+M/Z 


n;  m„  n,  ,  x 

nj°j  =  2~VieniVi - -  Ve{niOi) 


J 


mi  ne 


j*i 


f  m  ^ 

1/  J - Ly 

is  es 

s*e,i\  mi  J 


Ukk. 

mi  ne 


veYjnivi 

V  j*i 


which  may  be  rewritten 

».Z 


m„ 


v*+— v,s 
mi 


_0% 

ni 

n , 

-Hk 

m, 

(  ^ 
VeHnj»J 

J 


-Hl 


~  me 

Z  V;„ - -V. 


\ 


\  J 
\  ( 


m„ 


mi  j 
\ 


I n,v, 


Vis  +  ■  es 

s*e,i\  J 


V, 


Using  the  above  expression,  we  may  rewrite  Eq.  (2.5)  as 


Define 


1+i 

n„ 


\mi  vi 


V, 


= 


Vi 


k_ 

m; 


b-v(»,r.) 

n0 


{, 

l+Ok 

g-b  +  «/X( 

me  y 

vis+  —  ves 

ni 

vs+  — 

(  rr, 

V  V 

vie  ve 

l 

.  mi. 

s±ej\ 

mi  J 

ne 

l  mi 

h*i 

(2.9) 


/?  -  v.+^~ 


au 

\mi 


V -2v,„ 


me 
He  mi 


M  L 


n:  m0 

ij  '  ^  ej 

n0  m; 


r-z 


ni  me 

y.  +  — - - —  y 

r  in  1  K 

mi 


•/*' 


n.  me 

V  +— * - -  V”  . 

ij  v  ej 

ne  mi 


1-2— 

v  «,y 

/ 


n,  m 

Ke+——Ki 

ne  mi 


1-  — 
V 


(2.10) 


(2.11) 


where  we  have  once  more  made  use  of  the  relation  «,w,  v1(,  =  nemevei.  Noting  that  corresponding 
electron  and  ion  collision  frequencies  are  not  very  different  in  magnitude,  and  neglecting  mjmi 
compared  to  terms  of  order  unity,  we  may  approximate 
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2>,„+ive+(i -i 

n  j*i  V  "e 


(2.12) 


Note  that  for  a  single  ion  species  the  ion-electron  collision  term  vanishes  exactly.  In  the  absence 
of  a  background  gas  of  neutrals,  vf -»  0  and  the  concept  of  ambipolar  diffusion  becomes 
invalid.  This  may  be  a  concern  at  high  altitudes  where  neutral  densities  are  very  low  and  H+  is 
the  dominant  ion. 

With  the  definition  of  vf  the  expression  for  the  ion  flux  becomes 

nlo,=-^F\-—\b  V(nlT,)  +—b-V(n,Tl,)  +n,  1+^  g-b 

*1  “'l  J  L  m']  (2.13) 

,„A  mi  )  "A  '  h*‘  J 

Neglecting  me/mi  with  respect  to  terms  of  order  unity,  the  equation  simplifies  to 

■  V(»,7;)+ '  H”J.)  +  »,gb 

vf  [  mt\_  ne  J 

+n£Jvinvn+niYJvijVj+ni  vie  — — ve  L—»j  ' 

n  j*i  \  m<  e 

(The  relative  magnitudes  of  vie  and  ve  =  ven  +^y.vey  is  not  immediately  apparent,  so  both 
terms  are  retained.) 

Now  we  express  the  pressure  gradient  terms  in  terms  of  the  coordinates  Q  and  X  of  Anderson 
[1973].  First,  we  write  the  explicit  definition  of  b: 

b  =  —  ^Brr  +  B0Q  +  B^)  =  -smIr-cosIcosdQ-cosIsmd^ 


with  5  =  |B|  and  the  standard  definitions 


r  Br 
sin  7  =  — 

sin  8- — r== 
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Note  that  the  definition  of  i,  in  Anderson  [1973]  has  the  opposite  sign  so  that  i,  B  <  0.  With  our 
definition,  g  •  b  =  -gr  •  b  =  g  sin  /  where  g  =  |g| . 

With  the  above  definitions,  the  parallel  component  of  the  gradient  operator  is 

,  „  .  T  d  coslcosS  c  cos/sin<5  c 

bV  =  -sin/ - 

dr  r  30  rsin#  dtp 

Define  Q  =  yj (r0g,° )  where  y  is  the  scalar  potential  such  that  B  =  -  V?  .  Then 


b  V 


ccX 

3Q3X 


where  C  is 


and  X  is  defined  as 


.  TcO  cos /cos S  cQ 

C  =  sin/— + - —  + 

dr  r  30 


cos  1  sin  6  cQ 
rsin#  d(j) 


sinh(TQ) 

sinh(r^max) 


where  Qmax  is  the  value  of  Q  at  the  first  (northernmost)  point  on  the  field  line. 

For  GTIM,  we  assume  a  dipole  field  and  work  in  dipole  coordinates  (that  is,  in  a 
coordinate  system  aligned  with  the  dipole  axis).  In  that  case 


y{r,0,<p)  =  g,°  r0  f — 1  cos# 


Q(r,0,<p)  =  cos 6 

■  -2g°^j  cos#,  Be  =  g?^-j  sin#,  =0,  and  #  =  g,°^j  >/l+3cos2# 


so  that 


.  2  cos  6  . 

sin  /  =  —f=  =  and  sin  o  =  0 

Vl  +  3  cos2  6 
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In  any  case,  we  obtain  the  following  expression  for  the  ion  flux: 


m  = 


Hjg  sin  I  | 

V? 


k  1  3X  d_ 
+  m,.  vf  3Q[_3X 


+  vf  +  vf  +  vf  [V/e  m,  Ve ) %n 

C—l—faT; )+  ———{neTe ) 

3Q\_3Xy  ,J  ne3XXe  e’\ 


(2.14) 


Since  vf  contains  the  ratio  n,/we,  every  term  in  Eq.  (2.14)  is  non-linear  except  in  the 
special  case  of  a  single  ion  (n,  a  ne).  We  argue  that  non-linear  effects  produced  by  vf  are  small. 
If  «,  »  nM  the  last  term  of  Eq.  (2.12)  is  negligible  and 

vf  (»/»»;*) 


h  y*j 


On  the  other  hand,  if  «,  ~  write 


T  ”  nimi 


rij  neme 


nej  n,m, 


ZV  f  i  ^ 

V-  +  >  v,7  +  1 — L  — 
^  7  n  m 

n  j*i  \  e )  rni 


using  «,  ~  ny  ~  ne  and  m)  ~  mr  But  vw  ~  vy/  and  ^  « 1,  so  the  last  term  is  again  negligible. 
Finally,  if  ni  « itj  for  some  j,  then 


v'< 

R  yVi  V  J 

~'ZVi*+'LVJi  +  Vi* 


n  j±i 


and  the  dependence  on  «,/«e  disappears.  Thus,  in  all  cases,  ni/ne  has  only  a  small  effect  and  can 
simply  be  treated  as  a  constant  during  each  time  step. 


Now  consider  the  pressure  gradient  terms: 


PG.-i lwc% 

ntjVf  3Q  3c  n„  3c  n„TT;  3c' 


ne  3c 


(2.15) 


Define 
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(2.16) 


(,)  _  k  oX 
3  =  ~  m,vf  cQ 


so  that 


PG  s  -B{i] 


8c  ne  d K  HejTidt 


(2.17) 


To  cast  this  in  a  form  similar  to  that  used  in  LOWLAT  and  GTIM,  write  Eq.  (2. 1 7)  as 


PGs-EkH 


d_ 

8k 


(  n.  ' 

Ti+-^Te 

\  ne  J 


W, 


ne%8c 


d_ 
e  8c 


'n' 


\neJ 


(2.18) 


In  order  to  use  the  LOWLAT/GTIM  solver,  we  need  to  write  the  pressure  gradient  terms  as 


PG  m  -Bf 


^■(Tfn)*(F-G)n, 


(2.19) 


1  ~  (  \ 
where  T?f  =  T,  +  —  T  ,  F  =  —  Y  —  (n ,Te),  and  G=Te—  .  We  will  take  T? ,  F,  and  G  as 

n  n  ^  Asr  \  J  cr  *  Ar\  wi  I 


ne-£8c 


8c  [  n 


\ne) 


constant  during  a  time  step.  Let  us  examine  the  validity  of  this  approach. 


Case  1 :  ni  »  n /V/ 

In  this  case  ~  1  so  holding  and  G  constant  is  reasonable,  while  F  is  small. 

Case  2:  w,  «  ne 

In  this  case,  Tf®  »  7J,  and  G  is  very  small.  In  F,  ne  «  ny  for  some  /  *  /  so  neglecting  its 
variation  during  the  time  step  is  also  reasonable. 

Case  3:  n,  «  «y  (for  some  j) 

In  this  case,  the  approximation  is  not  valid,  but  we  can  hope  that  it  only  effects  a  few  gird 
points  in  the  transition  from  one  dominant  ion  to  another  (e.g.,  from  0+  dominant  to  H+ 
dominant). 

Now,  we  rewrite  the  ion  flux  as 
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(2.20) 


where 


4° = ^ +prZ  w. +-jrl  w ' 

vi  vi  n  yi  j*i 


+zsrvi>-—v’ 


™i  )j*i  n 


xrnJ 

±.—°j 


Now,  we  need  to  take  the  divergence  of  this  flux. 

V  •  (n,u,b)  =  b  •  V(n,uI)+  n ,-o,.V  •  b 


V  •  (n,-W/b)  =  +  ».-u.‘V  • b 


The  full  continuity  equation  becomes 


.cX  c 


%  =  *,  -  A,»,  +  C^(w)  -  »,«,V  •  b  -  „,V  •  V, 

^  OX 


where  V,  is  (E  x  B)/B2 .  Substituting  from  Eq.  (2.20) 


%  =  *i+K-U,.V-b-V.Vjn,+C.g|:  (2.24) 


(2.21) 


(2.22) 


(2.23) 


°r  f  „  -| 

f  -  *.  -1*1  *  v  •  VJ, + £(W] 

Keeping  in  mind  that  T=1  and  Gi  =  n„  we  can  make  the  following  identifications.  First  we  wnte 
the  multi-ion  diffusion  equation  in  the  notation  of  Anderson  [1973]: 


(2.25) 


+B(>7x  ^+Bllx  ^  ^ + B*)n‘ + Bg 


Then,  we  identify  the  coefficients  B\' 
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*,=i; 


52=-C^; 


(0-  k  d(')  ■ 


B^  = 


m,.vf  2 


4°  =  T^; 


#  -  -^+^rUl  +WZ;>>  r4(i|; 


2^°  =  ^0V-b; 


57=- 


dt  Q* 


4°  =  -K-  +  V  •  vx  +  -S5  V  •  b] ;  4°  =  a,; 


with  vf  =  v(>.  +%-  vie-^-v\  ve  =  X +  Z v*/ ’  U\\  =U-b  =  o„  and 

ne  mi  J  n  j 

„  .  rdO  cos  I  cos 6  dO  cos  I  sin  6  dQ  „7  ,  .  .  ,.  . 

C  —  sin  /  — — - + - — + - — .  We  have  assumed  that  all  neutral  species  move 

dr  r  dO  rsin#  d<f> 

together  with  a  common  bulk  velocity. 


Neutral  wind  models  give  the  wind  velocity  components  in  geographic  coordinates  rather 
than  geomagnetic  dipole  coordinates.  The  conversion  between  the  two  coordinate  systems  is  not 
described  in  this  document. 


To  avoid  ambiguity,  we  also  note  that  the  diffusion  equation  as  used  in  LOWLAT  and 
GTIM  is  written  in  a  slightly  different  form.  The  correspondence  between  Bk  as  used  above  and 
the  “COEF&”  that  appear  in  the  code  itself  is 

COEF 1  =  B\  COEF2  =  B2  COEF3  =  4°  COEF4  =  4° 

COEF5  =  4°  |COEF6  =  gj  |cOEF7  =  ~4'1  COEF8  =  4° 

COEF9  =  4° 


3.  Collision  terms. 

Schunk  and  Nagy  [1980]  provide  expressions  for  the  collision  terms  required  in  GTIM. 
We  have  adopted  these  expressions  except  of  the  resonant  charge  exchange  between  O  and  0+. 
For  this  process  we  have  adopted  the  formula  of  Pesnell  et  al.  [1993]. 

In  the  following  expressions  concentrations  (n)  are  in  m'3,  temperatures  are  in  K,  masses 
are  in  kg,  and  collision  frequencies  are  in  s'1. 

*st  ~ 

ms  +  mt 

a.  collisions  between  0+  and  other  species 
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31 
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=5-45*  l°~s#- 

e,He+ 


=  5.45x1 0“5-^C- 

Te% 
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